Motion of the hydrogen bond proton in cytosine 
and the transition between its normal and imino states 
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The potential energy surface of the H13 proton in base cytosine of the DNA molecules is calculated 
ab initio at the Gaussian98 MP2/6-311G(d,p) level. Two potential wells are found. One corresponds 
to the normal cytosine, while the other corresponds to its imino tautomer. The bindings of the proton 
in these wells are stable enough against the thermo-disturbance. The motions of the proton in these 
wells are oscillations around the nearest nitrogen atom like the pendula, and may move far away from 
the nitrogen atom to form the hydrogen bond with other bases. The estimated tunneling probability 
of the H13 proton from one well to another well shows that the life time of the proton staying in one 
of these wells is about 6xl0 2 yr. It is too long to let tautomers of cytosine be in thermodynamical 
equilibrium in a room temperature gas phase experiment. The biological significance of these result 
is discussed. 
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I. INTRODUCTION 

More than 40 years ago, Per-Olov Lowdin pj proposed 
a program for the research of mutation, aging, and can- 
cer formation on the basis of DNA changes due to the 
tunneling of the hydrogen bond proton through the bar- 
rier separating its normal and ill configurations. It was 
unable to realize this attractive program at that time, 
since both the theoretical chemical method and the com- 
putational ability were not enough developed. Having 
ab intio method with its fairly developed software pack- 
age and high performance computers in hand, one is able 
to try it now. Many works have been done along this 
direction. References are some earlier examples. 

However, we do not yet see the previous report on the 
ab initio calculation about the potential surface and the 
corresponding motion of the hydrogen bond proton. 

The four bases (A, T, C, G) in a DNA molecule 
usually stay in their normal configurations, but occa- 
sionally transfer to their unusual imino/enol configu- 
rations with small probabilities. It is the transitions 



A<-> A: 



T<-> T c 



imino' A ' ' A enob Gh~^G enol' ' v imino' 

The imino/enol tautomers will be paired in the way 

^imino ' ' ' ' ^imino' ^enol ' ' an< ^ ^' ' ' ^enol' 

It makes the possible changes 



and C^ C; 



C---G<->T---A 



(1) 



of the base pairs in DNA during its double replications. 



For an example, the transition C^ C 
G 



imino 



changes the 

pair C- • G to Cj mm0 G , and the pairing A- • • Cj mmo 
produces pairs A- • • Cj mmo and G- • -C in a first repli- 
cation. A further replication produces pairs T- • A, 
Cj mm0 • • -A ,C- • G , and C- • G. In the first case, a pair 
T- • A appears at the position originally occupied by the 
pair C- ■ G . This is a change of heredity code, therefore is 



serious. It may influence the development and reproduc- 
tion of a living system, cause aging, cancer, and mutation 
for examples. To understand and research these kinds of 
processes at the level of quantitative theoretical physics 
would be highly desirable. 
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FIG. 1: The bases cytosine and guanine and their pairing by 
hydrogen bonds (dashed lines) The attached number is for 
specifying the atom in a base 

Figure 1 and 2 show the structure formulae of base cy- 
tosine and its imino tautomer, as well as their different 
pairings. The difference between these two tautomers 
of cytosine lies mainly in the different positions of the 
proton H13. This difference in turn causes the different 
distribution of hydrogen bonds, and finally causes the 
different pairings. Therefore, to clarify the potential sur- 
face of proton H13 and its motion under the influence of 
this potential is the key for understanding the transition 
Cimino between cytosine tautomers and the change 
C- • • G -> T • • -A of the base pair in DNA. 

Section 2 is devoted to finding the potential surface 
of the proton H13 (the hydrogen bond proton) for base 
cytosine by ab initio calculation. In section 3, we prove 
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FIG. 2: The bases cytosine-imino and adenine and their pair- 
ing by hydrogen bonds 




that the proton H13 is bound in both, normal and imino, 
configurations, and it is ready for forming the right hy- 
drogen bond of pairing. In section 4 we calculate the 
tunneling probability of the proton H13 between these 
two configurations and discuss its biological significance. 
Section 5 contains the conclusion and a discussion of the 
reliability of the results. 



II. POTENTIAL ENERGY SURFACE OF THE 
PROTON H13 FOR CYTOSINE 

We have calculated the single point energy (SPE) of 
the proton H13 ab initio by MP2 method in Gaussian98 
package with the 6-311G(d,p) base set. For every given 
position of this proton, positions of other atoms in cy- 
tosine are optimized by the minimization of the energy 
of the electron system. It causes a small change of these 
positions in cytosine. Using a coordinate system fixed 
on the known structure formula of cytosine, the obtained 
map of the potential energy surface of the proton H13 is 
shown in figure 3. In this figure, the origin of the coor- 
dinate system is at the position of Nil in figure 1, the x 
axis is along its direction Nll-Cl, and the y axis is per- 
pendicular to x axis and points to the right. The bluer 
is the point in figure 3 the lower is the potential at that 
position. The number attached on the line shows the sin- 
gle point energy along the contour. We see clearly two 
potential wells. The left one corresponds to the cytosine, 
and the right one corresponds to its imino tautomer. It 
explains the structure of these two cytosine tautomers 
shown in figures 1 and 2. 



III. BINDING OF PROTON H13 IN THE BASE 
CYTOSINE AND ITS IMINO TAUTOMER 

Around the bottom of a well, we fit the calculated sin- 
gle point energies by a quadratic form of coordinates and 
find the minimum of SPE by differentiation. The mini- 
mum SPE of the left well is -393.9535268 a.u and that 
of the right well is —393.952654 a.u. Moving the origin of 
the coordinate system to the position of minimum SPE, 



FIG. 3: Contour diagram for the potential energy of the pro- 
ton H(13) in cytosine. The energy is in atom unit (a.u.), and 
the length is in the unit of pico-meter (pm). 



the terms containing first degree of coordinate disappear. 
Finding the principal axes of the coordinate system by 
the standard method, we may transform the quadratic 
form into a sum of squares of normal coordinate. For the 
left well, we take a coordinate system (£, rf) with the £ 
axis along the direction N11-H13 in cytosine of figure 1, 
and r\ axis perpendicular to it. The normal coordinates 
(X, Y) is related to (£, rj) by a rotation. 



0.94 
0.34 



-0.34 
0.94 



(2) 



The SPE as a function of normal coordinates (X, Y) is 



E = -393.9535 + 0.088937F 2 , 



(3) 



in which the single point energy E is in unit of a.u., and 
the unit of the length is angstrom. We see that E is inde- 
pendent of the coordinate X at the bottom of this well. 
The potential surface here looks like a valley along X 
direction. The transformation J5J shows it is a direction 
near the line N11-H13. Along this direction, the proton 
is easier to go far away from the atom Nil to form a hy- 
drogen bond with 015 in guanine. It explains the pairing 
shown in figure 1. The oscillation is along the Y direc- 
tion perpendicular to this direction. It is an oscillation 
around the atom Nil in cytosine, like a pendulum 'hang' 
at Nil. The frequency is 



12 x 0.088937 x 27.21 x 1.60217733 x 10" 19 / .„ . 
" = » / ~ : — / (2tt) 



1.67231 x 10- 27 x (10- 10 ) z 
= 3.427 x 10 13 (Hz) , 

and the zero point energy (ZPE) is 
\hv = 0.0026a.u . 



(4) 



(5) 
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Adding this ZPE with the SPE at the bottom of the 
well, we find that the energy of the proton H13 is 
—393.9509 a.u.. It is lower than the lowest maximum 
SPE -393.8895 a.u. around the well. The proton H13 is 
therefore bound in the normal cytosine, with a binding 
energy 6=0. 0614a. u.. 

For the right well in figure 3, we take the coordinate 
system (x , y), with the x' axis parallel the x axis in 
figure 3, while its origin is moved to the bottom of the 
well. The normal coordinates (X, Y) is related to (x', 
y') by a rotation 



0.30 
0.95 



-0.95 
0.30 



The SPE as a function of X and Y is 



E = -393.9527+ 0.12878F 2 



(6) 



(7) 



E is independent of the coordinate X again at the bottom 
of the well, and the potential surface looks like a valley 
along the X direction once more. The transformation (J5J) 
shows it is a direction almost perpendicular to the line 
Nll-Cl. Along this direction, the proton is easier to go 
far away from the atom N2 of cytosine-imino in figure 2 
to form a hydrogen bond with atom N2 in adenine. It 
explains the pairing shown in figure 2. The oscillation 
is almost along this line. It is an oscillation around the 
atom N2 in cytosine-imino, like a pendulum 'hang' at 
this atom N2. The frequency is 



2 x 0.12878 x 27.21 x 1.60217733 x 10" 19 



1.67231 x 10- 27 x (10- 10 )^ 



:4.12 x 10 i3 (Hz) 




and the ZPE is 



1 



-hv = 0.00313a.u . 



(9) 



The lowest maximum SPE around two wells are at the 
same point, it is a point between these two wells. The 
SPE —393.8895 a.u. at this point is called the energy of 
the transitional state. The sum of the ZPE and the bot- 
tom SPE is the energy —393.9496 a.u. of the proton H13 
in this well. It is also lower than the transitional state en- 
ergy. The proton H13 is bound too in the cytosine-imino. 
The binding energy here is 6'=0.0601 a.u. 

The thermo-energy in a living system may be esti- 
mated to be fcgT^O.OOla.u. at the room temperature 
T ~ 300K. We see that b and b' are much larger than 
the thermal energy. It implies that both cytosine and 
its imino tautomer are stable against the thermal distur- 
bance. 



IV. TUNNELING BETWEEN TWO 
TAUTOMERS AND ITS BIOLOGICAL 
SIGNIFICANCE 

According to quantum mechanics, there is still a pos- 
sibility of the proton tunneling from one well to another 
well. We first find the path of lowest potential connect- 
ing the centers of these two wells. It is described by the 
equation 

y(x) = 1.1417352 + 0.71244583.x + 0.3508047a; 2 

+ 0.04422011a; 3 + 0.06419875a; 4 - 0.56534255a; 5 
+ 0.25248558a; 6 (10) 

in the coordinate system used in figure 3, with the length 
unit being changed to angstrom. The potential energy in 
atom unit along the path is fitted by the formula 

V(x) = -393.908877+ 0.103691a; - 0.102525a; 2 

- 0.3063312a; 3 + 0.102250a; 4 + 0.355250a; 5 

- 0.195806a; 6 . (11) 

The potential barrier between two wells of the proton H13 
is shown in Figure 4. According to the WKB formula, 
the penetration probability of the proton through this 
barrier is 



W = exp • 



: exp ■ 



ft - . 
2 



\J 2nip [V(x) - E] dl 



2m p [V(x) - E] 



hj a \ 

:exp[-54.85] w 1.50929 x 10~ 24 



dx 
(12) 



Here E — —393.9496 a.u. is the total energy of the proton 
H13 in cytosine-imino, a and b are the x coordinates of 
the proton to be determined by the condition E = V(x) 
with b > a, and s is the lowest potential path terminated 
by points x = a and x = b. 

The penetration probability in unit time is 



P = vW = 3.427 x 10 13 x 1.50929 x 10~ 24 Hz 



5.168 x 10 



-H g -1 



(13) 



Therefore, the mean life time for a proton staying in one 
of these wells is 6xl0 2 years. It means the time needed 
for reaching the thermo-equilibrium is more than six hun- 
dred years. One could not see this equilibrium of gas type 
cytosine and its imino tautomer in his one hundred year 
life. We had better to assume that the large ratio between 
populations of cytosine and its imino tautomer in nature 
is not a result of thermo-equilibrium, but rather the natu- 
ral selection resulting in biological evolution. From 113|) 
we also see, in one's one hundred year life time, about 
sixteen percent of his cytosine bases may have sponta- 
neously transferred into some imino form. Some of them 
may transfer back into normal form of cytosine, while 
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FIG. 4: The potential barrier between two wells of the H13 
proton 

some others, according to Lowdin's conjecture, may in- 
duce aging, mutation or cancer in further replications of 
DNA. Of course, cytosine as well as other bases of DNA 
molecules are not isolated in living systems. Their tran- 
sitions to 'ill' tautomers must be influenced by their en- 
vironment. To clarify this influence, one has to do more 
complicated calculations for a base-environment system. 
A deeper research along this direction would be interest- 
ing. 



including the possible positions and directions of their 
hydrogen bond formation. It makes us be confident that 
this method is applicable and effective in the bio-physics 
research, and predict the transition probability between 
normal and imino states of cytosine. The result is rea- 
sonable. Its biological significance is discussed above and 
worthy to explore further. 



There are a large number of atoms in a bio-molcculc. 
The energy of the molecule is contributed by atoms in 
it and by their interactions. The energy change of the 
molecule in the physical, chemical or biological reactions 
is usually a tiny part of the total energy of the molecule. 
To study the behavior and the energy change of the bio- 
molecule in these reactions, we need very precise calcu- 
lations to keep many digits of the value be correct. For 
an example, as we have seen at the beginning of section 
3, to compare the single point energy of the proton at 
the centers of two wells, we need at least first 6 digits 
of the calculated values be precise. This is not easy in 
such a huge amount of numerical calculations. For the 
results reported here, the general tendency of the poten- 
tial surface may be reliable, but the difference of SPE for 
specially chosen pair of points may not be so, since the 
first 5 most reliable digits are subtracted. On the other 
hand, the tunneling probability may be reliable, since a 
part of statistical errors may be canceled by each other 
in an integration. Anyway, to develop the more precise 
programs with specified upper bound of possible errors 
would be very desirable in further researches along this 
direction. 



V. CONCLUSION AND DISCUSSION 

The ab initio quantum mechanical calculation repro- 
duces the structures of cytosine and its imino tautomer, 



This Work was supported by National Nature Science 
Foundation Committee of China with grant numbers 
10305001 and 10475008, and by the High Performance 
Computing Center of China (Beijing). 



[1] Per-Olov Lowdin, Rev. Mod. Phys. 35 (1963) 724. 

[2] J. W. Drake, E. F. Aliens, S. A. Forsberg, R. M. 

Preparata and E. O. Greening, Nature 221 (1969) 1128. 
[3] R. V. Wolfernen, J. Mol. Biol. 40 (1969) 307. 
[4] M. J. Nowak, K. Szczepaniak, A. Barski and D. Shrgar, 

Z. Natreforsch C 33 (1978) 876. 
[5] M. Raszka, and N. O. Kaplan, Proc. Nat. Acad. Sci. USA 

69 (1972) 2025. 
[6] R. R. Shoup, H. T. Miles and E. D. Becker, J Phys Chem. 

76 (1972) 64. 

[7] M. D. Topal, J. R. Fresco, Nature (London) 263 (1976) 
285. 

[8] M. D. Topal, J. R. Fresco, Nature (London) 263 (1976) 
289. 

[9] W. Saenger, Principles of Nucleic Acid Structure 



(Springer- Verlag, New York, 1984). 
[10] R. Rein and F. E. Harris, J. Chem. Phys. 43 (1965) 4415. 
[11] C. Colominas, F. J. Luque, M. Orozco, J. Am. Chem. 

Soc. 118 (1996) 6811. 
[12] C. N. R. Rao, Jerus. Symp. Quant. Chem. Biochem. 5 

(1973) 107. 

[13] M. Geller and B. Lesyng, Biochem. Biophys. Acta 417 
(1975) 407. 

[14] I. Kulakowska, M. Geller, B. Lesyng, K. Bolewska, and 
K. L. Wierzchowski, Biochem. Biophys. Acta 417, (1975) 
420. 

[15] K. G. Rao and C. N. R. Rao, J. Chem. Soc. Perkin Trans. 
2 (1973) 889. 



